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We study the mean-field phase diagram of the two-dimensional (2D) Hubbard model in the Lieb 
lattice allowing for spin and charge density waves. Previous studies of this diagram have shown that 
the mean-field magnetization surprisingly deviates from the value predicted by Lieb’s theorem [1] as 
the on-site repulsive Coulomb interaction (U) becomes smaller [2]. Here, we show that in order for 
Lieb’s theorem to be satisfied, a more complex mean-field approach should be followed in the case 
of bipartite lattices or other lattices whose unit cells contain more than two types of atoms. In the 
case of the Lieb lattice, we show that, by allowing the system to modulate the magnetization and 
charge density between sublattices, the difference in the absolute values of the magnetization of the 
sublattices, mLieb, at half-filling, saturates at the exact value 1/2 for any value of U, as predicted 
by Lieb. Additionally, Lieb’s relation, mLieb = 1/2, is verified approximately for large U, in the 
n £ [2/3, 4/3] range. This range includes not only the ferromagnetic region of the phase diagram of 
the Lieb lattice (see Ref. [2j, but also the adjacent spiral regions. In fact, in this lattice, below or 
at half-filling, mLieb is simply the filling of the quasi-flat bands in the mean-field energy dispersion 
both for large and small U. 


I. INTRODUCTION 


Despite intense research in the last few decades, the 
2D Hubbard model in the square lattice has remained an 
open theoretical problem in the field of the strong corre¬ 
lated systems [30]. Although it is known that at half¬ 
filling, the spin dynamics of the 2D Hubbard model is 
described by the Heisenberg antiferromagnetic exchange 
term [5], there is no consensus regarding the ground state 
magnetic phase diagram of this model. In fact, even at 
the mean-field (MF) level, depending on the magnetic 
phases allowed, different authors obtain different dia¬ 
grams for the square lattice [6]. The traditional order¬ 
ings are ferromagnetism, antiferromagnetism and para¬ 
magnetism Emm. Later, spin spiral phases, a general¬ 
ization of the previous three, were introduced [12) . The 
MF phase diagram became even more complex with the 
consideration of spatial phase separation QMHj- 


The Hubbard model in decorated 2D lattices has also 
been extensively studied, motivated by the search for 
metallic (fiat-band) ferromagnetism. These decorated 
lattices fall into three categories: the Lieb [I], Mielke 
m and Tasaki lattices m- All of these lattices share a 
common feature: the presence of flat bands in the energy 
dispersion relation. In the particular case of the Lieb’s 
lattices, the flat bands result from the topology of the lat¬ 
tice, while in the case of Mielke and Tasaki lattices, the 
flat bands reflect longer-range transfer integrals in the 
system. One of the most representative examples of dec¬ 
orated lattices is the Lieb lattice, which can be obtained 
from the 2D square lattice, for example, by inserting an 
extra atom between every two nearest-neighbours (see 
Fig. la). Each unit cell (shaded rectangle in Fig. la) has 
3 atoms, one of each kind: A, B and C, whose relative 
occupation is depicted in Fig. |lb[ in the limit of no in¬ 
teractions. Fig. [lc] shows the energy bands of the Lieb 
lattice in this limit. Examples of materials whose struc¬ 
ture resembles the Lieb lattice include La2_, E Sr. 2 ,Cu04 


and YBa2Cu307, two well-known high-T c superconduc¬ 
tors with weakly coupled Cu02 planes ?18, 19j. 

A theorem by Lieb p] states that, in the particular 
case of bipartite lattices (i.e., lattices with two sublat¬ 
tices, A and B, such that each site on sublattice A has 
its nearest neighbors on sublattice B, and vice versa), 
the ground state is ferromagnetic at half-filling (n = 1, 
or one electron per lattice site), as long as the number of 
atoms of each sublattice is different. However, for exam¬ 
ple in the case of the Lieb lattice (a line-centered square 
lattice this ground state should be identified with 

ferrimagnetism (2T|. In fact, although each sublattice is 
indeed ferromagnetic, there is antiferromagnetic ordering 
between every pair of nearest neighbours [2] (see Fig. la). 

The magnetic phase diagram of the Hubbard model 
in the Lieb lattice was recently studied by us [2]. We 
showed that the mean-field magnetization per unit cell 
at half-filling, mueb = (|ms| + |mc| — ItoaD/ 2, assum¬ 
ing the particle density (n) of the tight-binding limit (see 
Fig. lb) and the same magnetization (to) in the whole 
lattice, surprisingly deviates from the value predicted by 
Lieb’s theorem [T] as the on-site repulsive Coulomb inter¬ 
action ( U ) becomes smaller, although these two assump¬ 
tions are common in mean-field studies [2l [221 - 125] . Lieb’s 
theorem predicts that the magnetization per unit cell, 
mLieb, is 1/2 for any U at half-filling. Fig. [2] shows both 
the mean-field phase diagram of the Hubbard model in 
the Lieb lattice, and the value of mLieb, using the mean- 
field results from Ref. 01 For the results to agree with 
Lieb’s theorem, one should have mLieb = 1/2 for any U. 
Although in the strong coupling limit (U t), the mean- 
field result satisfies Lieb’s theorem, it is far from correct 
near the tight-binding limit (U = 0). 


In this manuscript, we study the magnetic phase di¬ 
agram of the Lieb lattice allowing for different average 
occupations (n^, ns, and nc) and magnetization am¬ 
plitudes (rriA, tub, and me) in each sublattice. We 
find that with these new considerations, Lieb’s relation, 
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FIG. 1. (a) The Lieb lattice is a line-centered square lattice, 
comprising three sublattices, A, B, and C and having one 
atom of each type in a unit cell. The circles represent atomic 
nuclei and the arrows represent spins. At half-filling, the Lieb 
lattice is ferromagnetic within each sublattice but antiferro¬ 
magnetic overall. Moreover, at half-filling, the total spin per 
unit cell is 1/2, as predicted by Lieb [T. (b) Plot of the tight- 
binding (U = 0) particle density of each sublattice of the Lieb 
lattice, A, B or C, as a function of the total particle density, 
(c) Plot of the tight-binding dispersion relation e(k x ,k y ) of 
the Lieb lattice. The flat band is made up entirely of B and 
C orbitals. 

mLieb = 1/2, is satisfied for any U at half-filling, and 
satisfied approximately for large U, in the n £ [2/3,4/3] 
range, which includes not only the ferromagnetic region 
of the phase diagram of the Lieb lattice in Fig. [2] but 
also the adjacent spiral regions (note that, away from 
half-filling, TOLieb is no longer the unit cell magnetization, 
but gives the difference in the absolute values of the mag¬ 
netization of the sublattices). An important point in this 
result is that finite mLieb reflects the existence of quasi¬ 
flat bands in the mean-field energy dispersion. These 
quasi-flat bands are present, not only for small 17, but 
also for large U. In fact, below or at half-filling, mLieb is 
approximately the filling of the flat bands both for large 
and small U. 

Our mean-field approach follows that of Bach and 
Poelchau [26] (see also Refs. H3l and l27jl . The formal¬ 
ism with further mathematical details can be found in 
Ref. J2SJ 

The organization of this paper is as follows. We begin 
by presenting some key results for the Hubbard model 
and mean-field. Secondly, we revisit the tight-binding 
limit of the Lieb lattice. We then proceed to adding elec¬ 
tronic interactions to the Hamiltonian and calculating 
its mean-field counterpart. Finally, we show our results, 
discuss their meaningfulness, compare them to analyti¬ 


cal calculations and take conclusions. In Appendix A, we 
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FIG. 2. (a) Mean-field magnetic phase diagram of the Hub¬ 
bard model in the Lieb lattice, obtained in Ref. [2] and (b) 
difference in the absolute values of the magnetization of the 
sublattices, mLieb = (|m s | + |mcj — |r?vi|)/2, of the Lieb lat¬ 
tice at half-filling, using the mean-field results from the same 
reference, where the on-site magnetization is assumed to be 
the same on all sublattices. 

briefly outline the derivation which leads to the results in 
section II, and in Appendix B, we explain in more detail 
our method of calculating saddle points, which we used 
to obtain the results in section V. 


II. MEAN-FIELD METHOD FOR THE 
HUBBARD MODEL 

In this section, we adapt a key result presented as The¬ 
orem 4.14 in Ref. [25, This derivation was first done by 
Lieb and his colaborators Ezj and simplified by Bach and 
Poelchau [26]; in Appendix A we show an adaptation of 
the derivation in Ref. [25, Alternative approaches can be 
found in Refs. [25] and (25} The result consists of the 
following abridged derivation. 

We begin with the Hubbard Hamiltonian, given by 

H = t c t,<7°y^ + U ^2 (!) 

(x,y),° x 

Here, t is the hopping parameter between nearest neigh¬ 
bours and cj, a (cy jCr ) is the creation (annihilation) oper¬ 
ator of an electron on site x (y) with spin a =j\ j,. The 
letters x and y denote lattice sites, (x, y) stands for near¬ 
est neighbors and U is the on-site repulsive. The total 
number of particles is N. 

With the intent of finding the mean-field Helmholtz 
free energy, Fhf, associated to this Hamiltonian (note 
that we use F because we work with fixed number of 
particles; if we worked with fixed chemical potential, we 
would use the grand canonical potential instead), we first 
replace the interaction term, with the Hartree 

and Fock terms, 
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Replacing the averages by the mean-field parameters m and n (see Appendix A for details), it follows that the 
mean-field Hamiltonian, Fl(m,n ), corresponding to the Hamiltonian in Eq. [T] is 


H(m,n)=t 4,<t C v,<y + U Y 

{x,y),tT X 




(3) 


where s x and fi x are the spin and electron density op¬ 
erators at the site x, and fh x and n x are the respective 
mean-field parameters. 

The Helmholtz free energy F(m,n ) is calculated from 
H(rh,n ) using the partition function, Z(rh,n), 

F(m,n) = —In Z(fh,ri) = — In ^Tr 

= -4 Tr ( ln (!+ e ~ 0h )) + - n2 ) > 

" X 

(4) 

where 

hxcrya' — t&cra' T ^x * @era’ ) &xyi (h) 

Here, o aa i is the vector of Pauli matrices. 

The important result is that in the Hubbard model, the 
minimum of the Helmholtz free energy, Fhf, corresponds 
to a saddle point of its mean-field counterpart, F(rh,n), 

Fhf = minmaxf(m, n). (6) 

m n 

See Appendix A for a detailed discussion of the above 
relation. 

Computing the partial derivatives of F(m,n) with re¬ 
spect to m and n and setting them equal to zero, we find 
the self-consistency relations n x = (n x ) and rh x = ( s x ). 
In the case of the Hubbard model, solving the Hartree- 
Fock equations self-consistently is actually equivalent to 
finding a saddle point of the mean-field energy, i 7 ’(m, n). 
If one fixes the particle density, n on all sites of the lat¬ 
tice, then the mean-field calculation is reduced to finding 
a minimum of the mean-field energy with respect to the 
spin density, m. Note that these results are for finite 
temperatures, but they are still valid in the limit T —> 0 
in the case of the Hubbard model, as shown in Ref. 1261 

References Il3l and [281 go on to apply Eq. [6] to the com¬ 
putation of a mean-field magnetic phase diagram of the 
Hubbard model in the square lattice, imposing two re¬ 
strictions that we do not adopt in this work. Namely, 
they restrict magnetic phases to ferromagnetic (F), an¬ 
tiferromagnetic (AF), and paramagnetic (P), and addi¬ 
tionally impose homogeneous particle density throughout 
the square lattice, in which case the extremization of the 
free energy given in Eq. [g] is reduced to a minimization 
problem. However, in the case of the Lieb lattice, charge 
modulation occurs even in the tight-binding limit. For 


a mean-field calculation to yield the correct result in the 
tight-binding limit, this charge modulation needs to be 
taken into account. To the extent of our knowledge, this 
work is the first application of the result in Eq. [6] which 
uses both n and m to extremize F(m , n). To accomplish 
this, we use the generalized HF theory, which turns the 
original minimization problem for the Helmholtz free en¬ 
ergy F into a saddle-point problem for the mean-field 
Helmholtz free energy F(m,ri) (see Appendices). 


III. THE LIEB LATTICE IN THE 
TIGHT-BINDING LIMIT 

The Lieb lattice is a square lattice, with a quarter of 
its atoms removed in a regular pattern. Introducing a 
different creation operator in each sublattice, , B\ and 
C\ the tight-binding term of the Hamiltonian of such a 
model, H t , is given by J3D3J 

t E E [(4 , y B x , v + 4./4 + H e.; m 

+ (4,i/-®a,W-l + 4,y4-l,y + H.C.)] . 

L x ( L y ) is the number of unit cells along the x (y) direc¬ 
tion. The hopping terms in the first line are intra-unit 
cell and the remaining are inter-unit cell. Its eigenvalues 
originate three energy bands, one of which is flat. The 
dispersion relation for periodic boundary conditions is 

£± = ±2f^/cos 2 y +C0S 2 y, (8) 

for the two non-flat energy bands, where k a = 2im a /L a 
with n a = 0,1, • • • , L a and a £ {s, y}. The flat band is 
L x x L y -fold degenerate with zero energy. These three 
energy bands are shown in Fig. [Ic] The three branches 
intersect at the point ( k x ,k y ) = (n,n). Expanding the 
dispersion relation in Eq. [8] around this momentum, we 
find the Dirac cones e 2 = t 2 (k x + fc 2 ). The flat band is 
built up from B- and C-type orbitals, while the lower and 
upper bands involve all three lattices A, B, and C. This 
lack of uniformity in the distribution of the sublattices in 
the energy bands justifies the difference in the occupation 
numbers of the sublattices presented in Fig. |lb| 
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IV. INTERACTIONS AND MEAN-FIELD 

In this section, we add interactions to the tight-binding 
Hamiltonian of the Lieb lattice, H t , and reduce the quar- 
tic dependance of the resulting Hamiltonian on the cre¬ 
ation and destruction operators to a quadratic one, using 
the mean-field approximation |12| . 

The key difference between our approach and previous 
approaches is that we allow sublattices A, B, and C to 
have a different average occupation number each, while 
keeping the total number of particles of the system, TV, 
fixed [on each point of the (n, U) phase diagram]. We be¬ 
gin by defining an average particle density on each sub¬ 
lattice, 

riA = n + 5 a 

n B = n + S B (9) 

nc = n + Sc, 

along with the number of particles on each sublattice. 
For sublattice A, this would be Na = vaL, where L = 
L x L y is the number of unit cells (which in turn equals 
the number of sites A). To keep the number of particles 
equal to N, we apply the restriction Na+N b +Nq = N, 
which is equivalent to 

Sa + S B + Sc = 0. (10) 


tj (n + S A 0 0 \ 

Hg = — 0 n + S B 0 1 , (14) 

^ \ 0 0 n + Sc ) 

and 

U (m a 0 0 \ 

H m = -w- 0 m B e iq y 0 . (15) 

1 \ 0 0 m c e iq * ) 

The matrix Hhf above is written in the basis 

{% B kl C k » A k+2,T B k + 2ql C k+2q .;}> where the VeCt0r ( 1 = 

( Qx,Qy ) defines the spin orientation in the system, as in 
the works by Dzierzawa 2?] and Singh [[5T] . In this paper, 
we assume that the spin spiral wavenumber q remains the 
same as in Ref. [2] even though we allow the system to 
have spin and charge modulation. H t (k) is the matrix 
that corresponds to the tight-binding term of the Hamil¬ 
tonian. All other terms have correspondence with the 
interaction terms of the Hamiltonian in Eq. [3] Namely, 

T( m A + m l + m c)> 

x 

T X) n x > ^r(( n + $a) 2 + {n + S B ) 2 + (n + Sc) 2 ), 

X 

f n x n x -A ^diag(n + Sa, n + S B , n + Sc) = H s , 

X 


Setting Sa = S B = Sc = 0, we would obtain a lattice 
with its particles evenly distributed, which is what hap¬ 
pens in the usual 2D square lattice: all sites have the 
same particle density. Aside from total particle number 
conservation (Eq. 10) and motivated by the symmetry of 


the lattice, we impose is that S B = Sc- This gives the im¬ 
portant relation S B = Sc = —Sa/2 => n B = |(3 n — ua), 
which leaves us with one unknown with respect to which 
F(rh,n ) needs to be maximized. We chose to maximize 
with respect to Sa- 

In our case, we work at T = 0, so that the mean- 
field free energy of the system can be found by summing 
the mean-field energies of the lowest levels that the N 
particles can occupy (this is the usual Fermi sea). On 
each point ( n,U ), the total energy of the system, 
is obtained by adding the lowest N eigenvalues of the 
mean-field Hamiltonian Uhf n 23 cm El, 


F/hf = 


H t (k) + H s 


//,, 


( 11 ) 


Hi H t (k + 2q) + H s J ’ 

and then adding the diagonal terms 

{m,A + m B + m c - (n + S A ) 2 - (n + S B ) 2 - (n + 8 c ) 2 )- 

( 12 ) 

The smaller matrices that compose the Hamiltonian Hhf 
are 

0 t(l + e iky ) t(l + e ik *) 

H t (k) = ( t(l +e~ ik y) 0 0 

t(l + e~ ik *) 0 0 

' ( 13 ) 


-f Jim ■ s x -A — ^■diag(m J 4, m B e lqy , ) = H m . 

X 

( 16 ) 

The extra imaginary coefficients involving components of 
q arise from coupling sites on unit cells other than the cell 
labelled as ( x,y ). From this point forward, we consider 
t = 1, so that U is given in units of t. 


V. RESULTS AND DISCUSSION 

Our results consist of the values of via, itib an d n A 
which, for each pair (n,U) £ [0,2] x [0,20], correspond 
to a saddle point of 2?hf, the mean-field energy of the 
Lieb lattice (see Appendix B for a more detailed ex¬ 
planation on how to find these saddle points). From 
these three quantities, we can calculate me = m B , and 
n B = nc = |( 3n—nA)■ We do not impose different occu¬ 
pation or magnetization, we simply let the system choose 
the values which lead to a saddle point of the mean-field 
energy. Before performing the calculations for the Lieb 
lattice, we tested the algorithm for the Hubbard model 
in a square lattice and found that the occupations of all 
four sublattices (A, B, C, and D) were the same, i. e., 
Sa = S B = Sc = S B = 0, while the system chose to have 
two different magnetizations, mA = m B and m B = me, 
reproducing our results in Ref. 23 

The results for the Lieb lattice are in Fig. [3] In the 
following subsections, we discuss each region of interest 
in more detail. In particular, we study the low U and 
high U regions separately, and finally the near-half-filling 
region [n « 1). 
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FIG. 3. Raw numerical results of our mean-field approach, namely for (a) 71a, (b) ub, (c) uia and (d) ms. The bold red lines 
highlight the U = 0 and U = 20 edges of the plots, and the bold blue line at the center of each plot allows a clearer visualization 
of each plot at half-filling (n = 1/2). 


A. Results near the tight-binding limit (U — > 0) 


As expected, for low U, the average occupations of 

re¬ 


tire sublattices, tia and ns = nc (Figs. 3a and 3b 


spectively), approach those of the tight-binding limit, de¬ 
scribed in section III and plotted in Fig. |lb| As for the 
magnetization amplitudes, sublattice A is paramagnetic 
(rriA = 0) for any n and sublattice B displays a behaviour 
similar to that of m in Ref. [2] for low U. In other words, 
ms is finite between n = 2/3 and n = 4/3 and zero oth¬ 
erwise. As explained below, this difference is due to the 
fact that the flat band only involves B and C orbitals. A 
first-order perturbation analysis follows, explaining this 
result. 


Without the U perturbation, the system becomes a 
Lieb lattice without electron-electron interactions and 
displays the electronic density in Fig. |lb| The energy 


dispersion relation comprises three bands, as in Fig. lc 
which can be doubly occupied with no additional energy 
cost (because of the absence of U). Keep in mind that the 
dispersive bands comprise A-, B-, and C-type orbitals, 
while the flat band comprises only B- and C-type orbitals. 
The introduction of the repulsive first-order perturbation 
modifies the total energy of the system by shifting the en¬ 
ergy bands and adding the diagonal terms of Eq. [12] The 
dispersive bands are shifted by + F|a and the flat band is 
shifted by — FTi. Additionally, the flat band (which used 
to have zero energy) splits into two bands, separated by 
an amount proportional to UrriB- 


At zero filling, the energy bands are in their original 
(U = 0) position, because all S and m are set to zero 
(having at least one finite m would lead to higher energy 
due to the term in Eq. 121. As we insert electrons in 


the system, they occupy the lowest states in the lower 
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dispersive band, with n A = 2 ub = 2 nc■ Due to the 
perturbation, this slowly causes the flat band to shift to 
lower energy and the dispersive bands to shift to higher 
energy. Note that before the filling n ss 2/3, the system 
is able to have the lowest energy by displaying paramag¬ 
netism (m A = ms = 0), because, to first order, the only 
dependance of the energy on to is in the diagonal term of 
Eq. |12| This dependance is maintained for any n, in the 
case of m A . At n ~ 2/3, the lower dispersive band is al¬ 
most full and, due to the flat bands having shifted by the 
small amount — electrons begin to occupy the flat 
bands, rather than the dispersive bands. From here on, 
up to n = 1, the magnetization amplitude of sublattice B 
(ms) increases, in order to separate the flat bands into 
two, and decrease the energy of the lower flat band, which 
is the one being filled. Due to the high density of states 
in the flat bands (which are related to B and C atoms 
only), newly added electrons choose to occupy sublattices 
B and C, until both flat bands are filled, which occurs at 
n ss 4/3. Between n = 1 and n = 4/3, the magneti¬ 
zation ms decreases again, because the lower flat band 
is full and electrons are now occupying the higher flat 
band, which has energy proportional to ms- For filling 
n £ [1,2], the behaviour is symmetrical to that of the 
n £ [0,1] region. 


B. Results in the strong coupling limit (U t) 


In this subsection, we discuss our results for high U, 
which we can assume to be nearly identical to those at 
U —> oo, for two reasons. Firstly, through inspection 
of the plots in Fig. [3] one readily realizes that the be¬ 
haviour at U = 20 is approximately the same as, say, 
U = 15, and therefore should not change with an in¬ 


Eu-¥ oo — (m-A + m B + m C ~ n A ~ n B ~ n C ) + -g 


where we have reintroduced the diagonal terms of Eq. [12] 
This expression can be simplified using the symmetries 
mentioned above: (i) 5a + 5 b + 5c = 0, (ii) ns = nc 
and (iii) ms = me- Using these three relations and 


UL 2 n i 3 2 9 2 n \ UL 

Eu^oo = m~ A + 2 m B - -n A - -n + 3 nn A ) + — 


crease in U. Moreover, if a certain value of U is enough 
to impose ferromagnetism (and therefore, for n < 1, the 
energy bands are singly occupied with spins in the same 
direction), higher values of U will have no additional ef¬ 
fect. Secondly, the analytical results for U —> oo which 
we present in the following paragraphs are in agreement 
with the numerical results in Fig. [3] for U = 20. 

One important first remark is that, as can be seen from 
Figs. [3a] and [3b] the value of n A or ns for U = 20 
along the line n £ [0,1], is the same as for U = 0 along 
the line n £ [0,2], but divided by two. This is because, 
at high U, the inequivalence of sublattices is imposed 
solely by the tight-binding terms of the Hamiltonian (see 
Eq. 11). Therefore, the behaviour of n A and hb is the 
same as in the tight-binding limit, albeit with all spins 
equally aligned. In fact, at U = 0, the tight-binding 
bands become doubly occupied without any additional 
energy cost, while for high U this cost is so high that all 
tight-binding states (thus, all sublattices) become singly 
occupied before double occupancies are created. 

To study the U —> oo limit from a perturbation theory 
point of view, we begin by setting t = 0 in the Hamilto¬ 
nian in Eq. ED and taking the result as the unperturbed 
Hamiltonian. Its eigenvalues are 


f {n A ± m A ), 


E{n B ± m B ), 


(17) 


f(n c ± m c )- 

These are six flat bands, with L states each. Positive 
to and negative to give the same set of eigenvalues, so 
let us assume positive to, with no loss of generality. At 
n < 1, electrons occupy the three lowest energy bands: 
U(? ia — tua), |(bb — ms) and w( n C — me), so that the 
total mean-field energy of the system is given by 


Y. { ~ tia - tua) + Y. {nn - m B ) + ^2{n C - me) , 
.n a n b n c 


(18) 


performing the summations up to some fixed N A = Ln A , 
the total energy for large U and n < 1 becomes 


n A {nA — m A ) + (3 n 


n A ) 



n A ) ~ m B 


(19) 


where the derivatives are taken with respect to m A , ms 


We find the ground state energy by taking VF/y^oc = 0, 
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and riA ■ The result is the self-consistency ua = n-A and 
the relations via = tla and ms = Ub- These two rela¬ 
tions hold true for U = 20, as can be realized by compar¬ 
ing ha with niA, and ng with ms in Fig. [3] Going back 
to Eu-t oo and replacing ua = tua and ub = nis, we find 
that the three bands in Eq. 17 become degenerate with 
zero energy. Note that setting ua = tua and ng = ms 
does not lead to a minimum of Eu^-oo, but to a saddle 
point, as expected. Again, Hartree-Fock mean-field the¬ 
ory is not about finding minima, but rather about find¬ 
ing self-consistency, as explained in section II. Finally, we 
remark that having ua = tua and ng = mg is a conse¬ 
quence of a ferromagnetic ground state with the bands 
being filled with only one spin direction. 

The next step in the perturbation analysis is to intro¬ 
duce the hopping terms of the Hamiltonian of Eq. H3 as 
the perturbation. This perturbation lifts the degeneracy 
of the three zero-energy bands. Up to first order, one of 
them, comprising B- and C-type orbitals, retains zero en¬ 
ergy, while each of the other two bands are pushed to pos¬ 
itive or negative energy, proportionally to t, and comprise 
orbitals of the three types, A, B, and C. This new energy 
band configuration mimics that of the tight-binding limit. 
This justifies the relative filling of the sublattices at high 
U (Figs. |3a| and 3b), and therefore their magnetization 
itia = tia and rag = ng. We now have a lower band 
involving A, B, and C orbitals, an intermediate flat band 
built up from B and C orbitals, and a higher energy band 
involving all three types of atoms. It follows that, as we 
begin inserting electrons in the system, sublattice A fills 
at two times the rate of the B/C sublattices. At n = 1/3, 
the lower band (with energy proportional to —t) is full 
and ua = 1/4 = 2 ub = 2 x 1/8. Between n = 1/3 
and n = 2/3, the flat band is filled up using only sublat¬ 
tices B and C, therefore ua does not change. Finally, for 
n £ [2/3,1], sublattices A, B, and C become gradually 
half-filled as the band with energy proportional to +t is 
filled. For filling n £ [1,2], the behaviour is symmetrical 
to that of the n £ [0,1] region. 


C. Results near half-filling (n ~ 1) 


At half-filling, the ground state is a q = (7r, n) phase. 
Lieb’s theorem [1] states that in bipartite systems, such 
as our Lieb lattice, the ground state at half-filling is fer- 
rimagnetic [21] . In our case, this ferrimagnetic ordering 
is ferromagnetic within each sublattice, and antiferro¬ 
magnetic between every two nearest-neighbour sites (see 
Fig. la). Although Lieb’s theorem does not provide in¬ 
formation on the spin per site, it does state that in our 
system the quantity \(2ms — itia) should be equal to 
1/2. This served as one of the motivations for this work, 
as previous studies of the Lieb lattice using mean-field 
[2] failed to yield the value 1/2 at low U. One of the 
reasons for this was that the magnetization amplitude 
was the same on all sites, A, B, and C, and the relative 
occupation of the three sublattices in the tight-binding 
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FIG. 4. (a) Plot of the difference in the absolute values of 
the magnetization of the sublattices, m Lieb = /(2 ms — ma), 
as a function of n and U, using the data from the plots in 
Figs. [3c|and [3d] The bold red lines highlight the (7 = 0 and 
U = 20 edges of the plot, and the bold blue line at the center 
allows for easier visualization of the behaviour of |(2ms — 
?ua) at half-filling. According to a theorem by Lieb [T], the 
value of /(2m_e — tua) at half-filling is 1/2, a value which the 
plot shows to have been achieved by our mean-field approach, 
(b) Phase diagram illustrating the relative behavior of the 
magnetization of sublattices A and B, as a function of n and 
U. 


limit was used for finite U. 

One can plot the function |(2 ms — tua), using our 
results for itia and ms, in Figs. [3c] and |3d| Such a 
plot can be found in Fig. [4a| At exactly n = 1 (the 
bold blue line at the center of the plot), our general¬ 
ized Hartree-Fock approach succeeds in yielding the re- 














suit |(2 ms — tua) = 5 , thus verifying Lieb’s theorem. 
This leads to the conclusion that our mean-field study, 
allowing for modulation of in and n in the Lieb lattice, 
produces more accurate results than imposing the same 
magnetization and electronic density in the whole lat¬ 
tice. Moreover, this plot indicates that Lieb’s theorem is 
verified approximately in the wider n £ [2/3,4/3] range. 

The phase diagram given by Fig. [ib| was constructed 
using the relative behavior of the magnetization of sub¬ 
lattices A and B. The dashed lines indicate the bound¬ 
aries where the behavior of the magnetization of each 
sublattice changes. The magnetization WLieb is exactly 
1/2 at the central (red) dashed line [n = 1). For small U , 
within the interval 2/3 < n < 4/3, only sublattice B has 
finite magnetization. For large U, in the n £ [2/3,4/3] 
range, Lieb’s theorem is verified approximately. Note 
that this range includes not only the ferromagnetic re¬ 
gion of the phase diagram of Fig. but surprisingly 
also the adjacent spiral regions. 

One other theorem, by Lieb and collaborators, refers to 
the particle density in the sublattices [32]. This theorem 
states that in the bipartite Hubbard model, there are 
no charge density modulations at half-filling. In other 
words, it means that at half-filling all sublattices are 
equally occupied and therefore half-filled, ua = rig = 
nc = 1- Plotting ha and ns as a function of U at fixed 
n = l, using our results for ha and ub (shown in Figs. 


the lines ua = 1 and ub = 1- As discussed in Section 
|V B| of this paper, the behaviour of ha and Ub at large 
U consists of two compacted copies of the behaviour at 
zero U. This holds true for any bipartite system. Conse¬ 
quently, for large U, not only are there not charge density 
modulations at half-filling, but also there are no charge 
modulations at n = 1/2 and n = 3/2. 

The contrasting behavior in Fig|4b|for large and small 
U can be interpreted as a consequence of the modifica¬ 
tions of the mean-field energy dispersion as U increases. 
More precisely, for small U, one has a single quasi-flat 
band (twice degenerate with energy e ~ 0 ), while for 
large U, one has two non-degenerate one-particle quasi¬ 
flat bands well separated in energy (e ~ 0 and e ~ U). 
Note that the flat bands for large U are present in the 
exact solution of the Hubbard model in the subspace 
of eigenstates associated with saturated ferromagnetism. 
As one can conclude from Fig. [4a] the difference in the 
absolute values of the magnetizations of the sublattices 
grows only when a flat band is being filled and in fact, 
we can write 


3a and p3b] respectively), we checked that we obtained 


WLieb = -( 2 tob - m A ) 


filling of flat bands, ( 20 ) 


for filling n < 1, both for large and small U. For n > 1, 
one has the reflected behavior of n < 1 . 


VI. CONCLUSIONS 


In summary, we have studied the Lieb lattice using 
a mean-field approach and allowing for charge and spin 
density modulation. Although theory about the cor¬ 
respondence between Hartree-Fock self-consistency and 
saddle points of the mean-field energy is relatively old 
(20 years old), to the extent of our knowledge, this is 
the first time it is in fact applied to a system where 
charge modulation is known to occur. We have found 
that, in the limits of low interaction (U 0 ) and very 
high interaction (U —> 00 ) results agree with what one 
would expect. Namely, the relative occupation of sublat¬ 
tices A and B of the bipartite Lieb lattice (where sublat¬ 
tice A comprises the atoms with four nearest neighbours, 
and B denotes the remaining atoms) in the tight-binding 
limit coincides with the results in the literature (for in¬ 
stance, in Ref. We have also found that the profile of 
the relative occupation of the sublattices in both strong¬ 
coupling and tight-binding are analogous, and one can be 
inferred if the other is known. The argument is relatively 
simple. On the one hand, in the U = 0 limit, the energy 
dispersion is governed by the tight-binding terms of the 
Hamiltonian only and the energy bands in the case of 
the Lieb lattice are as depicted in Fig. [Tc] On the other 
hand, in the limit U t, the energy bands are separated 
by a very high energy gap (of the order of [/), the lower 
bands corresponding to no double occupancies and the 
higher bands to double occupancies. The gradual filling 
of the system is done by singly filling all sites and only 
then jumping to the higher bands and doubly occupying 
all sites. Each one of these two filling regimes follows 
the relative occupation of sublattices that occurs in the 
tight-binding limit. 

At half-filling, our numerical mean-field results ver¬ 
ify two important exact results for the Hubbard model. 
Firstly, one theorem [1] states that in bipartite lattices 
with more B-type atoms than A-type atoms, the to¬ 
tal spin per unit cell at half-filling is equal to WLieb = 
(|R| — |A|)/ 2 , where |x| denotes the number of x-type 
atoms a unit cell. In the case of our bipartite Lieb lat¬ 
tice, we get (2 — l)/2 = 1/2. Our results are plotted in 
Fig- @ At exactly half-filling we obtained the expected 
value 1 / 2 . Additionally, for large U, in the n £ [2/3,4/3] 
range, we found that WLieb ~ 1/2. Interestingly, this 
range includes not only the ferromagnetic region of the 
phase diagram of Fig. Sb but also the adjacent spiral 
regions. Secondly, another theorem [32] states that in 
the bipartite Hubbard model, there are no charge den¬ 
sity modulations at half-filling, that is, at half-filling all 
sublattices are equally occupied and half-filled. Our nu¬ 
merically calculated relative occupations of sublattices A 
and B of a bipartite Lieb lattice, shown in Figs. [3a] and 


3b are in agreement with this theorem. In addition, we 


found that, for large U, not only are there no charge 
density modulations at half-filling, but also there are no 
charge modulations at n = 1/2 and n = 3/2. 

Away from half-filling, we found that, for large and 
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small U, the difference in the absolute values of the sub¬ 
lattice magnetizations (wLieb) grows or decreases only 
when a flat band is being filled and furthermore, for 
n < 1, m-Lieb is given approximately by the filling of the 
flat bands. Note that mLieb is the unit cell magnetization 
in the case of the q = (n, 7r), the ferrimagnetic phase of 
the Lieb lattice found at half-filling. 

We suggest that much of the above discussion is valid in 
the case of other bipartite lattices with flat bands in the 
energy dispersion. In the case of non-bipartite lattices 
with more than two types of atoms in the unit cell, the 
analysis is more complex, because, for instance, an an- 
tiferromagnet configuration between sublattices may not 
be commensurate with the unit cell. Bipartite lattices 
of various geometries can be realized by manipulating 
quantum dot arrays [33] or cold atoms in optical lattices 

SI- 


APPENDIX A: MIN-MAX THEOREM FOR THE 
HUBBARD MODEL 

In this appendix, we discuss how the mean-field 
method should be applied to the Hubbard model if parti¬ 
cle density is not fixed and in particular, we justify Eqj6 
Our approach follows the method presented in Ref. [26 
We start by presenting the usual mean-field approach and 
then we explain, using the method presented in Ref. [2B1 
why this approach is somewhat misleading. 


A. The usual mean-field method 


We begin with the Hubbard Hamiltonian, given by 

H = t 'y ( c x,o- c y,CT T U 'y ) (21) 

(x,y),a x 


We then replace the interaction term with the Hartree 
and Fock terms (see Eq. [2]). These terms are obtained 
by considering that each fermionic operator only devi¬ 
ates slightly from its mean value, so that in a product of 
operators, we can neglect quadratic terms in these devi¬ 
ations. 

The operators ^c x -|- and |C Xj j, can be identified 
with the particle density operator n x and the opera¬ 
tor. For consistence with the notation in reference [28], 
the spin operators in this appendix are 


4 c t 

4 c t 

cjc^ 


s+ = ^(s 31 + is y ) 
s~ = l(s x + is y ) 
C [ C i = s z . 


( 22 ) 


The vector s x is the spin density operator on site x. The 
peculiarity of this definition of s + and s~ is the factor 1/2 
which is often included in operators s x and s v instead. 
Replacing this in the interaction term of the Hamiltonian 
in Eq. [2l] gives 


Uhf = 


t y ) c x,trCy,a + U ^4 ^ ^( s x) {fl x ) ^ + 2 ( '(s x ) • 


(x,y),° 


(23) 


We now replace the averages by the mean-field parameters m x and n x . These parameters can be identified with the 
mean values of the magnetization and particle density, respectively, upon extremization of the mean-field free energy, 
i.e., when the self-consistency equations are satisfied. The Hartree-Fock Hamiltonian becomes 


U 


H{m, n) y ( h xaya ic\. IJ Cy^ a i + y ( 


- «*) . 


(24) 


where 


T'xcrya' — t^crcr' “1“ ^ ' ^crtr') $xy • 


(25) 


I 

Here, a aa i is the vector of Pauli matrices. The function F{rh,n) is calculated from H(rh,n) using the partition 

function, 


F(jh, n) = — jj In Z(m, n) = — ^ In (Tr (e 

= -jTr (In (1 + e~P h )) +^EK~ «*) . 


(26) 


At this stage, one usually finds the minimum free en¬ 


ergy, Fhf, by minimizing F(m,n ) with respect to the 
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mean-field parameters and this would lead to the usual 
self-consistency equations [HH] . Clearly, this works for 
fixed particle density, but in the previous expression 
we allow for variable particle density and the respective 
quadratic term has a negative coefficient. If one imposed 
a minimization with respect to the parameter n x , con¬ 
vergence would not be achieved in a numerical approach, 
unless one limits the possible values of n x to a certain 
interval, in which case the result of the numerical mini¬ 
mization would lie at the boundary of this interval. This 
reflects the fact that one should not minimize with re¬ 
spect to the parameter n x , but instead maximize, as we 
explain in the next subsection. 


Hamiltonian (the mean field approximation is associated 
with the state, not with the Hamiltonian). However, cal¬ 
culating min 7 F(y) going through all 7 is not pratical, 
so one introduces mean-field parameters. In the next 
paragraphs, we introduce these parameters, following a 
derivation in Ref. l26l For now, let us assume zero tem¬ 
perature. 

In the case of the Hubbard Hamiltonian, the Hartree- 
Fock energy functional, E( 7 ), can be written as (Ref. [261 
Lemma 2) 

E( 7 ) = Tr[T 7 ] + uY J [(n,) 2 - (4) 2 ], (28) 


B. A different perspective for the mean-field 
method 

In this subsection, we present the mean-field approach 
which should be applied when one takes into account the 
possibility of non-uniform particle density in a lattice. 
For a rigorous proof see, for instance, Refs. [26l - i28l 

Our goal is to know the thermal equilibrium state of 
the system which minimizes the free energy. These states 
are defined in terms of density matrices. Having an ex¬ 
act free energy would require having an exact partition 
function, which in turn would require an exact diagonal- 
ization of the Hubbard model. The Hartree-Fock method 
provides an approximation to the exact equilibrium state, 
in terms of many-body states of non-interacting parti¬ 
cles, replacing the quartic terms of the Hamiltonian by 
one-particle potentials (which are adjusted to provide the 
best possible approximation). The free energy obtained 
in the Hartree-Fock method provides an upper bound 
to the exact free energy (a consequence of the varia¬ 
tional theorem). Since the particles are independent, the 
Hartree-Fock state can be written as a one-particle den¬ 
sity matrix, 7 ^ = ( c\cj) [26]. The objective of the mean- 
field method is indeed to find the minimum free energy 
in the set of free energies associated with the possible 
states of N independent particles (or equivalently, asso¬ 
ciated with the possible one-particle density matrices), 

F hf = minF(y) = rnin^y) - -^(l)]- (27) 

1 1 P 

where /3 is the einverse temperature. In the previous 
expression, no mean field parameters are present and 
the free energy is determined for each 7 using the exact 


where x labels the lattice sites and T is the matrix whose 
elements t xy are the transition amplitudes of an electron 
to move from site x to site y or vice versa. The averages 
of the electronic and spin densities at site x are 


(n x ) = Tr 


I • 7 


(s x ) = Tr (I<g> <?) • 7 


(29) 


respectively. 

The first step to introduce the variational parameters 
is to use the simple fact that 

x 2 > 2 xy — y 2 for all x, y £ R n , (30) 

where equality holds if and only if x = y. This implies 

a ; 2 = max(2 xy — y 2 ), (31) 

v 

and therefore, we can write 

(fix) 2 = max {2 (h x )n x - n 2 x ) (32) 

n x 

and 

- (s x ) 2 = min {to 2 - 2 (s x ) ■ fh x \. (33) 

where n x and rh x are, at this point, an arbitrary constant 
and vector. Inserting this into Eq. [28j the energy for a 
given state 7 at zero temperature assumes the form 

E("f) = min max {i?(ri, to, 7 )} , (34) 

rh n 

where E(n, rh , 7 ) is given by 


js 

s 

II 

4 

M 

5 

H 

t=l> 

1 

Ql>. 

H 


. \ all sites 




+ u Y n D - 

all sites 


(35) 


The aim of the mean-field method is to obtain The crucial point of the mean-field method is the possi- 

min E( 7 ) = min min max {E(n, rh, 7 )} . (36) 

1 1 rh n 
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bility of exchanging the order of the extremization, which 
is implicitly done in standard mean-field. In fact, if par¬ 
ticle density is fixed, one has 

minE(y) = min min {E(m, 7)} = 7)} , 

7 7 m m 7 

(37) 

since the order of minimization is irrelevant, and one re¬ 
covers the usual mean-field picture, in which one has to 
minimize the energy with respect to the mean-field pa¬ 
rameter rh. 

On the other hand, if the particle density n x is allowed 
to actually depend on x, in Ref. [2fT| it was shown that 
indeed 


min min max {E(n, rh, 7)} = min max min {E(ra,m, 7)} . 

7 rh n ffl n 7 

( 38 ) 

Here, we present a simple physical interpretation for this 
result. 

First, one should note that Eq. [35] for given m and n 
corresponds to the energy of independent particles sub¬ 
jected to one-particle potentials which depend on the 
given rh and n. Looking at the left-hand side of Eq. [38j 
after one has minimized and maximized E(n,fh, 7) with 
respect to rh and n respectively, the remaining minimiza¬ 
tion (with respect to 7) will lead, at zero temperature, to 
a filled Fermi sea associated to these one-particle poten¬ 
tials. These one-particle potentials might, for instance, 
generate imbalance between the number of up and down 
spins and lead to finite magnetization. 


Second, Eqs. 32 and 33 associated with minimization 
and maximization with respect to rh and n, imply the 
self-consistency equations 


(h x ) = n x 
(s x ) = rh 3 


(39) 


and therefore, computing the left-hand side of Eq. [38] 
can be interpreted as the following instruction: ’’Among 
the set of one-particle density matrices that satisfy the 
self-consistency equations, find the one with minimum 
energy”. As we said above, this minimum energy cor¬ 
responds to a filled Fermi sea at zero temperature, so 
in the set of states corresponding to filled Fermi seas, 
there is one that satisfies the self-consistency equations. 
This helps us understand the right-hand side of Eq. [38] 
in the following way: ’’Among the set of all filled Fermi 
seas, find (the energy of) the state which satisfies the 
self-consistency equations”. 

This argument can be generalized for finite tempera¬ 
ture [SB], taking into account the entropy contribution. 
In this case, instead of filled Fermi seas, one has the one- 
particle density matrices that minimize the free energy 
in the case of independent particles. If instead of fixed 
number of particles, one imposes a fixed chemical poten¬ 
tial (in which case the grand-canonical potencial would 
replace the free energy), this one-particle density matrix 
would become the Fermi-Dirac distribution function. 


APPENDIX B: HOW TO EXTREMIZE E H f 


In this appendix, we explain in more detail the algo¬ 
rithm we used for finding saddle points. 

Our first approach to extremize Erf was perhaps the 
most intuitive non-brute force one. It consisted in max¬ 
imizing Erf with respect to 5a, 5b, and 5c using our 
results for q x , q y , rriA and tub of the Hubbard model in 
a square lattice [23]. The lattice size was also kept at 
100 x 100. Due to the two restrictions imposed (fixed 
number of particles in the system and equivalence of B 
and C sites), we are left with one parameter as maximizer 
of Erf- We used 5a- We note that the minimum of Erf 
with respect to 5a was found to be negative infinity (us¬ 
ing the results in Ref. [23] for niA and ms as starting 
points). Fixing q x , q y , iba and ms means that the mag¬ 
netic phases are kept the same as our previous ones, only 
the occupation numbers in the sublattices change. How¬ 
ever, finding saddle points by starting from a minimum 
with respect to one direction (the m direction) and then 
finding the maximum with respect to an orthogonal di¬ 
rection (the n direction) turned out to diverge on most 
points of the phase diagram. In fact, the suggestion that 
saddle points of a function can be found by starting at 
a random point in the function, minimizing the function 
with respect to the minimizer variables, and then using 
those new points to maximize the function with respect 
to the maximizer variables, is false. As a matter of fact, 
this statement remains false even in a more generalized 
case. One might think that by successively minimizing 
and maximizing the function, one would eventually reach 
a saddle point. This is also not necessarily true. When 
trying to use this method to extremize Erf, one finds 
that, in most points of the diagram, the value of Erf re¬ 
sulting from the last maximization and the value of Erf 
resulting from the last minimization differ by orders of 
magnitude comparable to those of E HF itself. Addition¬ 
ally, the values of rriA, tub, 5a and 5d which extremize 
Erf (or so one thought) do not converge on each suc¬ 
cessive iteration, they alternate between several possible 
results. 


In order to better understand why this method may 
fail to find saddle points, one can consider a much simpler 
example, the function f(x,y) = sin(x) cos (y). This func¬ 
tion is periodic, oscillates between -1 and 1, and has in¬ 
finitely many minima, maxima and saddle points (Fig. [ 5 ]). 
Let us now assume that we want to find the saddle points 
of /(x, y) using the method of minimizing and maximiz¬ 
ing several times alternating between the two, using x as 
minimizer and y as maximizer. If we start with a non¬ 
saddle point, the algorithm fails to find a saddle point and 
instead, after 2 ou 3 steps, alternates between absolute 
minima and maxima. Even if we start at a saddle point, 
successive iterations will move away from it and back to 
jumping between maxima and minima (see Fig. 5aI. No 
initial point will allow this algorithm to converge to a 
saddle point, no matter how close it is to it. This is pre¬ 
cisely what happened with the case of Erf- Depending 
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FIG. 5. (a,b) Contour plots of the function function f(x,y) = sin(a;) cos(y) in the range x, y G [—2tt,2-k\. The maxima are at 
the center of the red circles, the minima are at the center of the blue circles, and the saddle points are at the intersections of 
the green lines. The straight bold lines overlapped with the contour plots are the path taken by an algorithm which attempts 
to find a saddle point of f(x,y), by starting at point (—0.1,7r — 1.5) (red asterisk), and (a) successively minimizing in the x 
direction and maximizing in the y direction, and (b) successively minimizing in the direction which makes an angle of 7r/10 
with the x axis and maximizing in a direction which makes an angle of 7r/10 with the y axis, (c) Contour plot of the function 
f(x, y) = 2x 3 + 6 xy 2 — 3 y 3 — 150a: in the range (a;, y) G [—2,6] x [—2, 5]. The small (blue) dots represent minima of f(x, y) with 
respect to the x direction and the large (red) dot represents the maximum of these minima, which is a good approximation of 
a saddle point. 


on the pair (n, U) in question, £hf may or not behave 
in a way that allows saddle points to be found using this 
method. 

One possible alternative approach to finding saddle 
points of f(x,y ) is to attempt to extremize the func¬ 
tion in a different direction (or rotating the axes, which 
is another way to see it). Instead of using x and y as 
variables, we can use two auxiliary variables which are 
linear combinations of x and y but still orthogonal. For 
instance, we can alternate between the following two: 

• replace x —> Xq + r cos 6 , and y —l yo + r sin 9 , and 
minimize with respect to r; 


• replace x —> xq + r sin 9, and y —i yo — r cos 9 , and 
maximize with respect to r. 

Setting 9 = 0 would reduce this to simply minimizing 
with respect to x and maximizing with respect to y, as 
described above. Due to the symmetry of the simple 
function we used, we should ideally use 9 = 7 r /4 for the 
fastest convergence. For angles close enough to 7 r /4 we 
can indeed find a good approximation for a saddle point. 
Nevertheless, if we deviate too much from 7 r/ 4 , saddle 
points may not be found anymore. Using a ’’wrong” angle 
will make the algorithm diverge even if we start very 
close to a saddle point, even if the starting point is a 
saddle point itself. The algorithm will often be circling 
the saddle point (see Fig. 5b). 

The conclusion is that alternating between maximiza¬ 
tion and minimization is too sensitive to the initial point 
to be used to find saddle points. In the case of the sim¬ 
ple highly-symmetric function f{x,y), the saddle point 
is located at the center of the squares or rectangles that 
are drawn when connecting the points given by the algo¬ 
rithm, but this was shown to not be the case with Ahf- 


Moreover, this ideal angle depends on the pair (n, U ) 
which is being studied, and some angles may not even 
produce closed polygons, but rather diverge to infinity in 
a certain direction. 


The working alternative is failproof in the sense that 
it will always converge to one saddle point. It consists of 
calculating the value of f(x, y) for many points in an area 
which is known to contain at least one saddle point, mak¬ 
ing a list of the maximum value of the function for each 
x and then finding the minimum of all the maxima that 
were found. If f[x 1 y) is continuous we will end up on a 
saddle point. In the case of Uhf, we have not two vari¬ 
ables but three: rriA, tub and 5a- In order to achieve an 
acceptable precision, it would be reasonable to calculate, 
say, 100 values of E^f in each direction (tua, tub , and 
5a), for a total of 10 6 values for each pair (n,U) in our 
phase diagram, which could take a long time and would 
have a precision of about two decimal places. A more 
efficient alternative is to divide each direction into fewer 
parts, say 7 or 8 , finding an initial approximation to the 
saddle point, and then repeating this a dozen times con¬ 
sidering a smaller hypercube centered on the new point. 
Assuming each direction is divided into 7 parts (i.e. we 
calculate 8 function values in each direction) and the side 
of the hypercube is halved with each of 12 iterations, we 
now calculate a total of 7 3 x 12 sa 4000 values of U H f per 
pair (n, U). This was the procedure we used, which has a 
precision of around 1/4096, or three decimal places. An 
illustration of this method is shown in Fig. [5cJ where we 
have used the function f(x, y) = 2x 3 + 6 xy 2 — 3 y 3 — 150x 
and successfully found the saddle point (3,4). 
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